Radiation-MHD models of elephant trunks and 
globules in H II regions 
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Abstract: We study the formation and evolution of pillars of dense gas, known as elephant trunks, at the 
boundaries of H II regions, formed by shadowing of ionising radiation by dense clumps. The effects of mag- 
netic fields on this process are investigated using 3D radiation-magnetohydrodynamics simulations. For a 
simulation in which an initiaify uniform magnetic fiefd of strength |B| ~ 50 fiG is oriented perpendicufar to the 
radiation propagation direction, the field is swept into alignment with the pillar during its dynamical evolution, 
in agreement with observations of the "Pillars of Creation" in M16, and of some cometary globules. This effect 
is significantly enhanced when the simulation is re-run with a weaker field of ~ 18 /iG. A stronger field with 
|B| ~ 160 ;uG is sufficient to prevent this evolution completely, also significantly affecting the photoionisation 



process. Using a larger simulation domain it is seen that the pillar formation models studied in iMackey & Lim 



(120101) ultimately evolve to cometary structures in the absence of dense gas further from the star. 
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1 Introduction 



The interstellar magnetic field in the Eagle Nebula (Ml 6) was measured by ISugitani et al.l (l2007h us- 
ing near infrared polarisation observations of background stars. They found that the field orientation 
within the massive pillars of gas and dust in Ml 6 dHester et al.l . ll996t) . known as the "Pillars of Cre- 
ation", is aligned with the long axis of the pillars and with the UV radiation propagation direction, 
but misaligned with the ambient magnetic field (see Fig. [3] below). If the pillars have formed dynam- 



^ t^ig.lJjDeiow). it tne pmars nave lormea aynam- 
Williams. Ward-Thompson. & WhitworthI boOlt 



ically due to shadowin g of ionising radiation (e.g. ^ . _ ^ , _ _ 

Lim & Mellema . 2003 ) then it appears the field has been re-oriented by this process. ISugitani et al.l 
(120071) suggest that this may constrain the ambient field strength in M16 because ionised gas pres- 
sure appears to do minate the dynamics. Here we test this suggestion and build on our earlier non- 
magnetised results (|Mackey & LimL 12010 *. hereafter ML 10) by using 3D radiation-magnetohydrodynamics 
(R-MHD) simulations to model the formation of pillars and globules due to shadowing of ionising 
radiation by pre-existing dense clumps. 

A uniform grid, 2nd order accurate ^ finite volume code (se e ML 10) is used for the MHD calcula- 
tions with a Roe-type Riemann solver (|Cargo & GalliceLll997|) . The short characteristics ray-tracer is 
used to track direct monochromatic ionising radiation and the on-the-spot approximation for diffuse 
radiation. Radiative cooling is calculated explicitly for recombining hydrogen and using a cooling 
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Table 1: Simulation parameters for model R5, showing grid resolution (zones), domain coordinates 
in parsecs (Xmin, X^ax) relative to the radiation source which has monochromatic photon luminosity 
and photon energy as indicated. Clump positions, peak number density (tih), Gaussian scale radius 
(ro) and total mass (M) are as indicated; the background magnetic field (in /zG), number density (nn) 
and gas pressure (pg) are similarly indicated. 



Name 


X 


y 


z 


Further information. 


Zones 


384 


256 


256 


Uniform Cartesian grid. 


-^min 


1.5 


-1.5 


-1.5 


position in parsecs relative to the source. 


-^max 


6.0 


1.5 


1.5 


position in parsecs relative to the source. 


Source 











L^ = 2x 10^°s~\ huo - 13.6 eV = 5.0 eV 


Clump 1 


2.30 








nil = 10^ cm-^, ro = 0.09 pc, M = 28.4 M0 


Clump 2 


2.75 





0.12 


nil = 10^ ro = 0.09 pc, M = 28.4 M0 


Clump 3 


3.20 





-0.12 


nn = 10^ cm-^, ro = 0.09 pc, M = 28.4 Mq 


Background 


Bo = [0, 


0, 


53] 


nn = 200cm^^,pg = 1.38 x 10^^^ dynecm"^ 



curve for other elements, with neutral/molecular gas cooling treated very approximately by exponen- 
tial cooling with cooling time tc = 10 kyr (see MLIO, cooling mo del C2). This code has been exten- 
sively testecfl and checked for consistency \yith previous work (e.g. Falle et al. 1998 ; Lim & Mellema 
2OO3I : iMellema et aDl2006l : iKmmholz et~aDl2007l : [Hennev et al.ll2009l MLIO and references therein). 



2 Simulations 



We ha ve performed a number of 3D R-MHD simulations, to be presented in more detail in lMackey & Lim 



(|2011h . Here we present results from one of the simulations, denoted R5. This model consists of a 
monochromatic point source of ionising photons placed 2.3 pc from the nearest of three dense gas 
clumps (Gaussian density profiles) in a uniform background medium (riu = 200 cm^^). Simulation 
parameters are described in Table [1] the same model was presented on a smaller simulation domain 
with zero field in MLIO but here a uniform magnetic field of B = [0, 0, 53] fiG is imposed in the 
initial conditions. 

The evolution of dense gas over 400 kyr is shown in Fig. [H where the initial field is horizontal 
and in the plane of the images. The evolution is broadly similar to that of the R-HD model in MLIO 
for the first ~ 250 kyr, but there are notable differences. Radiation-driven implosion (t ^ 100 kyr) is 



more pronounced along field lines than across them, leading to a flattened structure (cf. i Henney et al. 



20091) : the narrow tail in the second panel is actually a sheet seen edge-on and the clumps are also 
flattened. Following this (200 — 300 kyr) magnetic pressure generates a stronger re-expansion and 
hence the pillar/globule is larger and has a lower density than the R-HD case. The re-expansion leads 
to fragmentation in this model (t > 400 kyr), just beginning in the right-most panel. It is seen that the 
simulation ultimately evolves from being pillar-like (dense elongated structure) to a cometary globule 
(dense flattened head, low density neutral tail) in the absence of dense gas further from the star. We 
have confirmed that this is also true when the model is run with zero magnetic field (i.e. identical to 
model 17 in MLIO but with a larger simulation domain). 



For more information, including test results, see |http : / / www .astro . uni-bonn . de/ ~ jmackey/ jmac/[ 
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Figure 1 : Volume-rendered images of gas density in the 3D R-MHD simulation R5 described above, 
shown at times t = 25, 100, 200, 300, 400 kyr. Tick-marks at 0.25 pc intervals show the physical 
scale (numbers refer to the tick-marks immediately below them). The logarithmic number density 
scale (cm^^) is shown above, with the opacity function used for image generation as the grey curve 
(gas with riH ^ 10^ cm^^ is transparent). Only part of the simulation domain is shown in each panel. 
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Figure 2: Column density (colour scale in cm~^) and projected magnetic field (lines show orientation 
only) for the R5 simulation at 250 kyr (above) and 400 kyr (below). The projection LOS is n = 
X sin(20°) +y cos(20°) in the simulation coordinate system, so the initial magnetic field (B = 532 /iG) 
is vertical and ± to the LOS. The horizontal axis labelled P is therefore P = x cos(20°) — y sin(20°) 
i.e. a 20° rotation from the x-axis about the 2;-axis. 
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3: Figure 9b from ISugitani et al. 
showing near-IR absorption po- 
larimetry observations of the magnetic field 
orientation in the central pillar in Ml 6 (© 
the Astronomical Society of Japan; used 
with permission). North is to the top, East 
to the left, and the image width is ^ 2'.5. 
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Figure 4: Ratio of volume-averaged parallel field 
{\Bx\) to perpendicular field {^jBy + Bf) as a 
function of time in simulations R2, R5, and R8. 
Red line: average of all cells (R5); Black: average 
of only dense cells with nn > 5000 cm~^ (R5); 
Green: simulated with |B| 3x weaker (R2); Blue: 
with |B| 3x stronger (R8). 



3 Discussion 



The projected magnetic field orientation was calculated by a density-weighte d integration of "Sto kes 
parameters" Q and U for the perpendicular magnetic field along the LOS (cf. [Arthur et all |2010|) : 
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with TT-max = 2.5 X 10 cm~'^. Here x — y h the image plane and z is distance along the LOS 
(i.e. observer rather than simulation coordinates). The projected field is then recovered from {Q) 
and ([/) using trigonometric relations. Formally this has units of /xG cm"^, but the normalisation is 
irrelevant for the orientation. The projected field orientation is over-plotted on column density maps 
in Fig. 121 projected such that the initial field is fully in the image plane and vertical. Ait = 250 kyr the 
clumps have undergone significant dynamical evolution and the pillar-like structure is slowly being 
flattened due to the rocket effect. The magnetic field, while relatively unaffected in ionised gas and 
the low density tail region, has clearly been dragged into alignment with the pillar in dense gas due 
to the pillar's acceleration away from the star. At 400 kyr only a flattened cometary globule remains, 
but the field is even more strongly aligned with the acceleration direction. This clearly shows that an 
initially perpendicular field is swept into alignment with the pillar during its dynamical evolution, and 
that this alignment remains once the structure evolves to a cometary morphology. This is in agreement 
with the observations o f the pillars i n M 16 (reproduced here in Fig. [3l ISugitani et al.ll2007h and some 
cometary globules (e.g. iBhatt et al.Ll2004|) which also have the field aligned with the cometary tail. 

The behaviour of the same simulation with a magnetic field 3 x weaker (R2) and 3x stronger 
(R8) is rather different however (to be discussed in more detail in iMackey & Liml l201l|) . For R2 
with S ~ 18 //G the field is more easily deformed by photoionisation-induced gas motions and the 
alignment is even stronger, whereas for R8 with B ~ 160 //G the field is dynamically dominant 
and barely changes from its initial perpendicular state. This is shown in the right panel of Fig. |3] 
where the ratio of the volume averaged field parallel and perpendicular to the radiation propagation 



direction (x) are plotted as a function of time for R2, R5, and R8. The evolution of dense gas with 
TT-H > 5000 cm"^ is shown for all three models, and for the full simulation volume only for R5 (R2 
and R8 are very similar). It is seen that while the ratio does not evolve significantly when averaged 
over the full simulation volume, in dense gas the situation is very different. For R5 the parallel field 
increases rapidly during the implosion phase and there is a subsequent slower increase during the 
acceleration phase to roughly equal field strength in both components. This trend is stronger for 
the weak field model R2, whereas for the strong field model R8 the field orientation remains almost 
constant throughout the si mulation. This differen ce between the weaker and stronger field simulations 
is larger than was found by iHenney et al.l (|2009|) and deserves further investigation. 



4 Conclusions 



The results presented here show that both radiation-driven implosion and acceleration of clumps by 
the rocket effect tend to align th e magnetic field with the radiation propagation direction in dense 
neutral gas. As was suggested by ISugitani et al.l (|2007b . the effectiveness of this alignment is depen- 
dent on the initial field strength. For simulations which have similar gas densities and pressures to 
conditions in Ml 6, we have shown that a magnetic field of strength ~ 160 fiG is sufficient to prevent 
any significant field reorientation. These results suggest that an ambient field of |B| < 160 fiG (and 
more likely |B| ~ 50 fiG) is required to explain the observed field configuration in the Ml 6 pillars if 
the pillars formed via the mechanism we ar e modelling. Detai l ed R-H D models have recently been 



performed with dynamic initial conditions (iGritschneder et al.L |2010|) rather than the initially static 
models of ML 10, allowing a potentially more realistic comparison with observations. Addition of 
magnetic fields to their simulations would be very useful to assess how significantly a non-uniform 
initial magnetic field will impact on the results presented here. 
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